t1 <- Sys.time()

sampd <- list.dirs("/media/knausb/E737-9B48/releases/scaffolded/")
sampd <- grep("scaffolded/$", sampd, value = TRUE, invert = TRUE)
sampd <- grep("EDTAOutput", sampd, value = TRUE, invert = TRUE)

#sampd <- sampd[ c(2:5, grep("EH23", sampd)) ]
#sampd <- sampd[2:5]

#sub(".+//", "", sampd)
sampn <- unlist(lapply(strsplit(sampd, split = "//"), function(x){x[2]}))

Queries

my_queries <- data.frame( 
  query = c("45S", "26S", "5.8S", "18S", "5S", 
            "CS-237_satelite", "CsatSD_centromere_237bp", "CsatSD_subtelomere_370bp"),
  Length = c(9382, 1544, 157, 1552, 259, 557, 237, 370 )
  )
knitr::kable(my_queries)

query	Length
45S	9382
26S	1544
5.8S	157
18S	1552
5S	259
CS-237_satelite	557
CsatSD_centromere_237bp	237
CsatSD_subtelomere_370bp	370

nuccomp

#i <- 1

# for( i in 1:length(sampn)){
#   cat( "Sample: ", sampn[i], "; ", i, " of ", length(sampn), ".\n", sep = "" )
#   cmd <- paste("~/gits/fasta2nuccomp/fasta2nuccomp.py ",
#                sampd[i] ,
#                "/",
#                sampn[i],
# #               "/",
#                ".softmasked.fasta.gz",
#                sep = "")
# 
#   if( !file.exists( paste(sampn[i], ".softmasked_nuccomp.csv", sep = "") ) ){
#     system( cmd )
#   }
# }

get_nucs <- function( sampd ){
  cat("__Processing chromosomes.\n")
  sampn <- unlist(lapply(strsplit(sampd, split = "//"), function(x){x[length(x)]}))
  cmd <- paste("~/gits/fasta2nuccomp/fasta2nuccomp.py ", 
               sampd,
               "/",
               sampn,
#               "/",
               ".softmasked.fasta.gz",
               sep = "")

  if( !file.exists( paste(sampn, ".softmasked_nuccomp.csv", sep = "") ) ){
    cat( "Running nuccomp: ", sampn, ".\n" )
    system( cmd )
  }
  
  nucs <- read.csv( paste( sampn, ".softmasked_nuccomp.csv", sep = "") )
  nucs <- nucs[ grep( paste(sampn, ".chr", sep = ""), nucs$Id ), ]

  if( length( grep("chrY", nucs$Id) ) == 0 ){
    last_chrom <- nrow(nucs)
    nucs <- rbind(nucs, NA)
    nucs[ last_chrom + 1, -1] <- 0
    nucs[ last_chrom + 1, 1] <- "chrY"
    nucs[ last_chrom + 1, 2] <- 1
  }
  if( length( grep("chrX", nucs$Id)) == 0 ){
    last_chrom <- nrow(nucs)
    nucs <- rbind(nucs, nucs[last_chrom, ])
    nucs[ last_chrom + 0, -1] <- 0
    nucs[ last_chrom + 0, 1] <- "chrX"
    nucs[ last_chrom + 0, 2] <- 1
  }
  return( nucs )
}

fast_win

#i <- 1

# for( i in 1:length(sampn)){
#   cat( "Sample: ", sampn[i], "; ", i, " of ", length(sampn), ".\n", sep = "" )
#   cmd <- paste("~/gits/hempy/bin/fast_win.py ",
#                sampd[i] ,
#                "/",
#                sampn[i],
# #               "/",
#                ".softmasked.fasta.gz",
#                " --win_size 1000000",
#                sep = "")
# 
#   if( !file.exists( paste(sampn[i], ".softmasked_wins.csv", sep = "") ) ){
#     system( cmd )
#   }
# }

get_wins <- function( sampd ){
  sampn <- unlist(lapply(strsplit(sampd, split = "//"), function(x){x[length(x)]}))
  
  cat("__Processing windows.\n")
  
  cmd <- paste("~/gits/hempy/bin/fast_win.py ", 
               sampd ,
               "/",
               sampn,
#               "/",
               ".softmasked.fasta.gz",
               " --win_size 1000000",
               sep = "")

  if( !file.exists( paste(sampn, ".softmasked_wins.csv", sep = "") ) ){
    cat( "Running nuccomp: ", sampn, ".\n" )
    system( cmd )
  }
  
  wins <- read.csv( paste( sampn, ".softmasked_wins.csv", sep = "") )
  wins <- wins[ grep( paste(sampn, ".chr", sep = ""), wins$Id ), ]
  wins$chrn <- sub(".+chr", "", wins$Id)
  wins$chrn[ wins$chrn == "X" ] <- 10
  wins$chrn[ wins$chrn == "Y" ] <- 11
  wins$chrn <- as.numeric(wins$chrn)
  
  # Scale and center.
  # Scaline of CG windows (chromosome width) is on a per chromosome basis.
#  wins$CGs <- 0
  wins$CGs <- wins$CG / wins$Win_length * 100
#  wins$notCG <- (wins$Win_length / 1) - wins$CG
#  wins$notCGs <- 0
  for( j in unique(wins$Id) ){
    wins$CGs[ wins$Id == j] <- wins$CGs[ wins$Id == j] - min(wins$CGs[ wins$Id == j], na.rm = TRUE)
    wins$CGs[ wins$Id == j] <- wins$CGs[ wins$Id == j] / max(wins$CGs[ wins$Id == j], na.rm = TRUE)
    #
    # wins$CGs[ wins$Id == j] <- wins$CG[ wins$Id == j] - min(wins$CG[ wins$Id == j], na.rm = TRUE)
    # wins$CGs[ wins$Id == j] <- wins$CGs[ wins$Id == j]/max(wins$CGs[ wins$Id == j], na.rm = TRUE)
    #
    # wins$notCGs[ wins$Id == j] <- wins$notCG[ wins$Id == j] - min(wins$notCG[ wins$Id == j], na.rm = TRUE)
    # wins$notCGs[ wins$Id == j] <- wins$notCGs[ wins$Id == j]/max(wins$notCGs[ wins$Id == j], na.rm = TRUE)
  }
  wins$iCGs <- 1 - wins$CGs

  #wins$ATs <- 1 - wins$CGs
  wins$ATs <- wins$Win_length/1e6 - wins$CGs

  genes <- read.table( 
    paste(sampd, "/", sampn, ".primary_high_confidence.gff3.gz", sep = "" ), 
    sep = "\t" )
  genes <- genes[genes[, 3] == "gene", ]

  # Windowize
  wins$gcnt <- 0
  for(i in 1:nrow(wins)){
    tmp <- genes[genes$V1 == wins$Id[i] & genes$V4 >= wins$Start[i] & genes$V5 < wins$End[i], ]
    wins$gcnt[i] <- nrow(tmp)
  }

  # Scaling of gene count windows is on a per assembly basis.
  wins$gcntsc <- wins$gcnt - min(wins$gcnt)
  wins$gcntsc <- wins$gcntsc / max(wins$gcntsc)
  my_index <- round( wins$gcntsc * 100 )
  my_index[ my_index <= 0] <- 1

  # Color ramp.
  #wins$gcol <- heat.colors(n=100)[ my_index ]
  #wins$gcol <- colorRampPalette(c("yellow", "orange", "red"))( 100 )[ my_index ]
  #wins$gcol <- colorRampPalette(c("red", "orange", "yellow"))( 100 )[ my_index ]
  #wins$gcol <- colorRampPalette(c("#0000FF", "#228B22", "#A0522D"))( 100 )[ my_index ]
  #wins$gcol <- colorRampPalette(c("#87CEEB", "#3CB371", "#228B22", "#A0522D"))( 100 )[ my_index ]

  #wins$gcol <- viridisLite::plasma(n = 100, alpha = 1, begin = 0.1, end = 1)[ my_index ]
  wins$gcol <- viridisLite::magma(n = 100, alpha = 1, begin = 0.2, end = 1.00)[ my_index ]
  
  return(wins)
}

Plot

#i <- 1
#i <- 2

#class(nucs$Length)
#nucs$Length

# wins <- read.csv( paste( sampn[i], ".softmasked_wins.csv", sep = "") )
# wins <- wins[ grep( paste(sampn[i], ".chr", sep = ""), wins$Id ), ]
# 
# # Scale and center.
# wins$CGs <- 0
# for( j in unique(wins$Id) ){
#   wins$CGs[ wins$Id == j] <- wins$CG[ wins$Id == j] - min(wins$CG[ wins$Id == j], na.rm = TRUE)
#   wins$CGs[ wins$Id == j] <- wins$CGs[ wins$Id == j]/max(wins$CGs[ wins$Id == j], na.rm = TRUE)
# }
# wins$chrn <- sub(".+chr", "", wins$Id)
# wins$chrn[ wins$chrn == "X" ] <- 10
# wins$chrn[ wins$chrn == "Y" ] <- 11
# wins$chrn <- as.numeric(wins$chrn)
# #wins$ATs <- 1 - wins$CGs
# wins$ATs <- wins$Win_length/1e6 - wins$CGs
# #wins[1:3, ]
# #tail(wins)
# # cwins <- wins[wins$CGs == 1 & !is.na(wins$CGs) , ]
# 
# genes <- read.table( 
#   paste(sampd[i], "/", sampn[i], ".primary_high_confidence.gff3.gz", sep = "" ), 
#   sep = "\t" )
# genes <- genes[genes[, 3] == "gene", ]
# 
# # Windowize
# wins$gcnt <- 0
# for(i in 1:nrow(wins)){
#    tmp <- genes[genes$V1 == wins$Id[i] & genes$V4 >= wins$Start[i] & genes$V5 < wins$End[i], ]
#    wins$gcnt[i] <- nrow(tmp)
# }
# 
# wins$gcntsc <- wins$gcnt - min(wins$gcnt)
# wins$gcntsc <- wins$gcntsc / max(wins$gcntsc)
# my_index <- round( wins$gcntsc * 100 )
# my_index[ my_index <= 0] <- 1
# 
# 
# # Color ramp.
# #wins$gcol <- heat.colors(n=100)[ my_index ]
# #wins$gcol <- colorRampPalette(c("yellow", "orange", "red"))( 100 )[ my_index ]
# #wins$gcol <- colorRampPalette(c("red", "orange", "yellow"))( 100 )[ my_index ]
# #wins$gcol <- colorRampPalette(c("#0000FF", "#228B22", "#A0522D"))( 100 )[ my_index ]
# #wins$gcol <- colorRampPalette(c("#87CEEB", "#3CB371", "#228B22", "#A0522D"))( 100 )[ my_index ]
# 
# #wins$gcol <- viridisLite::plasma(n = 100, alpha = 1, begin = 0.1, end = 1)[ my_index ]
# wins$gcol <- viridisLite::magma(n = 100, alpha = 1, begin = 0.2, end = 1.00)[ my_index ]
# 
# #wins$chr <- sub("^.+chr", "chr", wins$Id)
# #
# #wins[1:3, ]

Ideo function

#plot_ideo <- function(nucs, wins) {
#sampd <- sampd[1]

plot_ideo <- function( sampd ) {
  suppressPackageStartupMessages(require(ggplot2))

  sampn <- unlist(lapply(strsplit(sampd, split = "//"), function(x){x[length(x)]}))
  cat( "Processing ", sampn, ".\n" )
  nucs <- get_nucs( sampd )
  wins <- get_wins( sampd )
  
  blstn <- read.csv( 
#    paste("/media/knausb/Vining_lab/knausb/blast_projects/todd_rrna/blastn/", sampn, "_blastn.csv", sep = ""),
    paste("/media/knausb/Vining_lab/knausb/blast_projects/todd_rrna/blastn_uc/", sampn, "_blastn.csv", sep = ""),
    header = FALSE)
  colnames(blstn) <- c("qseqid","qlen","sseqid","slen","qstart","qend",
                       "sstart","send","evalue","bitscore","score","length",
                       "pident","nident","mismatch","positive","gapopen",
                       "gaps","ppos","sstrand") 
  #
  #length(grep("chr", blstn$sseqid, invert = TRUE))
  
  blstn <- blstn[grep("chr", blstn$sseqid), ]
  blstn$chrn <- sub(".*chr", "", blstn$sseqid)
  blstn$chrn[blstn$chrn == "X"] <- 10
  blstn$chrn[blstn$chrn == "Y"] <- 11
  blstn$chrn <- as.numeric(blstn$chrn)  
  
  # blstn[1:3, ]
  # table(blstn$qseqid)
  
  chr_df <- data.frame(
    start = 1,
    end = nucs$Length,
    chr = nucs$Id
  )
#  chr_df$chrn <- sub(".+chr", "", nucs$Id)
  chr_df$chrn <- sub(".*chr", "", nucs$Id)
  chr_df$chrn[chr_df$chrn == "X"] <- 10
  chr_df$chrn[chr_df$chrn == "Y"] <- 11
  chr_df$chrn <- as.numeric(chr_df$chrn)
  
  chrom_wid <- 0.02
  p <- ggplot2::ggplot()
  p <- p + ggplot2::geom_rect( 
    data = chr_df, 
    ggplot2::aes( xmin = chrn - chrom_wid,
                  xmax = chrn + chrom_wid,
                  ymin = start, ymax = end),
    fill = "#DCDCDC",
    color = "#000000"
  )

  thinw <- 0.28
  p <- p + ggplot2::geom_rect( 
    data = wins, 
    ggplot2::aes(
#      xmin = chrn - ATs * thinw,
#      xmax = chrn + ATs * thinw,
      xmin = chrn - iCGs * thinw,
      xmax = chrn + iCGs * thinw,
      # xmin = chrn - notCGs * thinw,
      # xmax = chrn + notCGs * thinw,
      ymin = Start, 
      ymax = End),
    fill = wins$gcol,
    color = NA
  )

  # Theme
  p <- p + ggplot2::scale_x_continuous( 
    breaks = chr_df$chrn,
    labels = chr_df$chr
  )

  p <- p + scale_y_continuous(
    breaks = seq( 0, 120e6, by = 10e6), 
    labels = seq( 0, 120, by = 10)
  )
  
  p <- p + ggplot2::theme_bw() + 
    ggplot2::theme( panel.grid.minor.x = ggplot2::element_blank(), 
           axis.text.x = element_text(angle = 60, hjust = 1),
           axis.title.x=element_blank(),
           panel.grid.major.y = ggplot2::element_line( linewidth = 0.4, color = "#C0C0C0", linetype = 1 ),
           panel.grid.minor.y = ggplot2::element_line( linewidth = 0.4, color = "#C0C0C0", linetype = 3 )
          )

  p <- p + ggplot2::ylab("Position (Mbp)")
  p <- p + ggtitle( sampn )

  my_locus <- ""
  table(blstn$qseqid)
  grep("5S_", blstn$qseqid)
  
  
  # Cannabinoids
  my_rows <- grep("AB2|LY|NC_044376.1:c427", blstn$qseqid)
  blstw <- 0.48
  p <- p + ggplot2::geom_rect( 
    data = blstn[my_rows, ], 
    ggplot2::aes( 
      xmin = chrn - blstw,
      xmax = chrn + blstw,
      ymin = sstart, 
      ymax = send),
    fill = "#228B22",
    color = '#228B22'
  )
  
  # 5S
  my_rows <- grep("5S_", blstn$qseqid)
  blstw <- 0.48
  p <- p + ggplot2::geom_rect( 
    data = blstn[my_rows, ], 
    ggplot2::aes( 
      xmin = chrn - blstw,
      xmax = chrn + blstw,
      ymin = sstart, 
      ymax = send),
    fill = "#000000",
    color = '#000000'
  )
  
  # 45S
  my_rows <- grep("45s_|26s_|5.8s_|18s_", blstn$qseqid)
  blstw <- 0.48
  p <- p + ggplot2::geom_rect( 
    data = blstn[my_rows, ], 
    ggplot2::aes( 
      xmin = chrn - blstw,
      xmax = chrn + blstw,
      ymin = sstart, 
      ymax = send),
    fill = "#B22222",
    linewidth = 2,
#    linewidth = 4,
    color = '#B22222'
  )
  
  # Cent, subteol
  my_rows <- grep("AH3Ma_CS-237_satelite|CsatSD_centromere_237bp|CsatSD_centromere_370bp", blstn$qseqid)
#  blstw <- 0.48
  p <- p + ggplot2::geom_rect( 
    data = blstn[my_rows, ], 
    ggplot2::aes( 
      xmin = chrn - 0.5,
      xmax = chrn - 0.2,
      ymin = sstart, 
      ymax = send),
    fill = "#1E90FF",
    color = '#1E90FF'
  )
  
  p
}

#p1 <- plot_ideo(nucs, wins)
#p1 <- plot_ideo(nucs, wins)
#p1

t4 <- Sys.time()

sampd <- list.dirs("/media/knausb/E737-9B48/releases/scaffolded/")
sampd <- grep("scaffolded/$", sampd, value = TRUE, invert = TRUE)
sampd <- grep("EDTAOutput", sampd, value = TRUE, invert = TRUE)

#sampd <- sampd[ c(2:5, grep("EH23", sampd)) ]
#sampd <- sampd[ c(2:11, grep("EH23", sampd)) ]

#sampd <- sampd[ c(1:2, grep("EH23", sampd)) ]
#sampd <- sampd[ c(1:10, grep("EH23", sampd)) ]

#sampd <- sampd[ c(1:16, grep("EH23", sampd)) ]
#sampd <- sampd[ c(1:22, grep("EH23", sampd)) ]
#sampd <- sampd[ 1:30 ]
#sampd <- sampd[2:5]

#undebug(plot_ideo)
# p1 <- plot_ideo( sampd[1] )
# p1 <- plot_ideo( sampd[5] )
# p1 <- plot_ideo( sampd[6] )
# p1

# plist[[1]] <- plot_ideo( sampd[1] )
# plist[[2]] <- plot_ideo( sampd[2] )
# 
# plist[[1]]
# plist[[2]]


plist <- vector( mode = "list", length = length(sampd) )
for( i in 1:length(sampd) ){
  message("Processing sample ", i, " of ", length(sampd), ".")
  plist[[i]] <- plot_ideo( sampd[i] )
}

# plist[[1]]

#sub(".+//", "", sampd)
#sampn <- unlist(lapply(strsplit(sampd, split = "//"), function(x){x[2]}))

library(ggpubr)

my_ps <- 1:4
#ggarrange( plotlist = plist[my_ps], ncol = 2, nrow = 3, labels = c("A", "B", "C", "D", "E", "F"))


for(i in 1:ceiling(length(plist)/4) ){
  my_ps <- 1:4
  my_ps <- my_ps + (i - 1) * 4
  if( my_ps[length(my_ps)] > length(plist) ){
    my_ps <- my_ps[1]:length(plist)
  }
  print(ggarrange( plotlist = plist[my_ps], ncol = 2, nrow = 2,
                   labels = c("A", "B", "C", "D")))
}

Figure X. Ideogram of the Cannabis sativa genome. Each chromosome is presented as a set of 1 Mbp windows where the width is proportional to the abundance of the ‘CpG’ motif (wide is low abundance, narrow is high abundance) and colored by gene density (cool or blue colors are low gene density and hot or yellow colors are high gene density). The 45S, 26S, 18S, and 5.8S probes are in firebrick red, the 5S probes are in black, Cannabinoids (CBCAS, CBDAS, THCAS, OAC) are in forest green, the centromeric and telomeric probes are in Dodger blue and are offset from the center of each chromosome.

# for(i in 1:(length(plist)/6) ){
#   my_ps <- 1:6
#   my_ps <- my_ps + (i - 1) * 6
#   if( my_ps[length(my_ps)] > length(plist) ){
#     my_ps <- my_ps[1]:length(plist)
#   }
#   print(ggarrange( plotlist = plist[my_ps], ncol = 2, nrow = 3, labels = c("A", "B", "C", "D", "E", "F")))
# }



#ggarrange( plotlist = plist[1:6], ncol = 2, nrow = 3, labels = c("A", "B", "C", "D", "E", "F"))
#ggarrange( plotlist = plist[7:12], ncol = 2, nrow = 3, labels = c("A", "B", "C", "D", "E", "F"))

t99 <- Sys.time()
t99 - t1

## Time difference of 5.075594 mins

Ideos

2023-02-26

Queries

nuccomp

fast_win

Plot

Ideo function